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We develop a self-consistent version of perturbation theory in Liouville space which seeks to com- 
bine the advantages of master equation approaches in quantum transport with the nonperturbative 
features that a self-consistent treatment brings. We describe how counting fields may be included 
in a self-consistent manner in this formalism such that the full counting statistics can be calculated. 
NonMarkovian effects are also incorporated. Several different self-consistent approximations are 
introduced and we discuss their relative strengths with a simple example. 
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I. INTRODUCTION 

Master equation (ME) techniques play an important 
role in quantum transport — not only in describing state- 
of-the-art experiments e.g. [iHEl, but also in contributing 
to the theoretical development of the field; topics such as 
full counting statistics (FCS), both at zero [ol and finite 
Q frequency; quantum coherent phenomena [3]; and non- 
Markovian (NM) effects have all been investigated 

with master equations. The central object in such ap- 
proaches is the Liouvillian or kernel, which generates the 
evolution of the reduced system density matrix in time 
and includes not only the internal dynamics of the sys- 
tem but also the effects of the coupling to the leads and 
to any further environmental degrees of freedom. ME Li- 
ouvillians can be derived in a number of different ways, 
from the standard textbook approaches, e.g. [H!, [I!], to 
more transport-specific ones such as the wavefunction ap - 
proach of Gurvitz [l3|, and the real-time diagrams [ij- 
\u\ . A reformulation of this latte r appro ach, the Liou- 
villian perturbation theory (LPT) [18l - l20 | . will form the 
basis of the current work. 

Despite their many advantages, standard ME ap- 
proaches are restricted to the regime in which system 
and environment are but weakly coupled, with the Li- 
ouvillian obtained as a perturbation series in the corre- 
sponding coupling strength. There is much interest in 
going beyond this situation and to develop nonpertur- 
bative approaches for the description of quantum trans- 
port beyond the ME. Some recent examples include equa- 
tions of motio n f2ll , renormalisation group [2^ and self- 
consistent [1^ |24| techniques. 

In this paper we develop a self- consistent (SC) version 
of the LPT, which seeks to combine the advantages of 
MEs with a nonperturbative treatment of the system- 
lead coupling such that effects such as, for example, 
level-broadening and higher-order tunnel processes are 
included. The method is self-consistent in the sense that 
system self-energies are calculated not with free propa- 
gators, but rather with ones that include the self-energy 
itself. We describe how counting fields may be included in 
a self-consistent manner in this formalism, which means 
that we are not restricted to the calculation of stationary 
current only (as in, e.g. [23j). but can, in principle, access 



the complete FCS of the system. NonMarkovian (NM) 
effects are also included in our theory from the outset. 

Within this framework, we introduce several differ- 
ent SC approximations. The simplest of these is just 
a SC version of the familiar sequential-tunneling kernel. 
This approximation is equivalent to the SC Born approx- 
imation f25| . appearing in a setting similar to that of 



Ref. [23[. By considering the single-resonant level (SRL) 
model as a test case, we see that although this type of 
kernel can provide an excellent description of the sta- 
tionary current even for large couplings (F/fesT 3> 1) it 
gives inaccurate results for the shotnoise, and is unsuit- 
able for the calculation of counting statistics in general. 
The problems of this kernel are remedied to some extent 
by the other two kernels described here, which contain 
not only sequential-type diagrams, but also cotunneling 
ones. With these kernels excellent results are obtainable 
for both current and shotnoise beyond weak coupling. 

This work is based on the exposition of the LPT given 
in Ref. H^, hereafter referred to as L To save needless 
repetition, the reader will be this paper for many of the 
details of the original LPT and the notation employed 
here. 



II. TRANSPORT MODEL 



We begin by specifying our general transport model 
with the Hamiltonian H — H^cs + -ffs + ^ composed of 
reservoir, system, and interaction parts. In its diago- 
nal basis the system part reads Ha, = '^^Ea\a){a\ with 
a) a many-body state of Na electrons. We assume a 
noninteracting reservoir Hamiltonian iJrcs = X]fe Q('^fca + 

f-a)(^\a^ka, with lead index a that includes spin and any 
other relevant quantum numbers; uJka is the energy of 
the fcth mode in lead a, aka a lead annihilation operator, 
and with /^^ the chemical potential of lead a. In equi- 
librium, the reservoir electrons are distributed according 
to the Fermi function /(w) = 1/ (e"/'=«'^ l), which, 
since we assume a uniform temperature, is the same for 
all reservoirs. Coupling between system and reservoirs is 
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described by the single-particle tunneling Hamiltonian 

V = ^ ^ tkaviO-\a^"^ ^" ^kam^ln'^ka, (1) 
ham 

where dm is the annihilation operator for single-particle 
level m in the system, and tkam is a tunneling amplitude. 



III. SELF-CONSISTENT LIOUVILLIAN 
PERTURBATION THEORY 

For completeness we give a brief review of the essential 
elements of the standard Liouvillian perturbation theory 
— for more details, the reader is referred I as well as 
to [isl . [l^ . Our starting point is to write down the von 
Neumann equation for the evolution of the total system- 
plus-reservoirs density matrix: 

p{t) = -nH,p(t)]=Cp{t). (2) 

This defines the Liouvillian super-operator C ~ 
—i[H,*]^ which, in accordance with the decomposi- 
tion of the Hamiltonian, consists of three parts: C = 

Crcs + Cs + Cv with Crcs = [-^res, • ], Cs = -i [Hs, • ], 

and Cv = —i[V,» ]. In terms of the Laplace-transformed 
total density matrix p{z) = dte~^* p{t), the solution 
of Eq. © is 

Piz) = J^pW- (3) 

Assuming a factorizable initial density matrix p(0) with 
reservoirs in equilibrium and system in state ps(0), trac- 
ing over reservoir degrees of freedom results in the fol- 
lowing reduced density matrix of the system: 

psiz) = Tr,.e. {p{z)} = _L_ps(0), (4) 

with the effective system Liouvillian W(z) = £3 + '^{z)- 
The self-energy T,{z) arises from the coupling with the 
leads and will be kept in its z-dependent nonMarkovian 
form. 

Equation is exact if kernel >V(z) is. In practice, 
however, the self-energy contains an infinite sum of terms 
and must be approximated in some way. In the straight- 
forward perturbative approach, the memory kernel is cal- 
culated as the series S(z) = E„^^"H-2): with n the num- 
ber of tunnel vertices in the term, and then approximated 
by simple truncation at a certain value of n. In previous 
works, expansion has been considered up to fourth order, 
such that the self-energy is approximated as 

with the first term describing sequential tunneling and 
the second, cotunneling. 

In describing the constituent terms of the self-energy 
and the various modifications to follow, it is useful to 



employ a diagrammatic representation. In this represen- 
tation, we write the sequential self-energy term as 

- G~~G. (6) 

z 

The translation of this diagram into an analytical expres- 
sion is discussed in Appendix lAl Here it suffices to note 
that the Gs represent tunnel vertices (of which there are 
two at sequential order) and these are connected by a sin- 
gle central line that corresponds to the free propagation 
of the system (i.e. under the action of free system Liou- 
villian >Cg). The overline connecting the tunnel vertices 
denotes a bath contraction between them. Such diagrams 
are similar to those of Ref . [l^ , but simplified such that 
we leave the indices implicit here. 

At fourth order in the perturbative expansion of the 
self-energy there are two diagrams, each with four tunnel 
vertices and three free propagators, but with differing 
patterns of contraction. In our diagrammatic notation, 
we write 



I I 1 I I 1 1 I 

y}^^^{z)^G-G-G-G^G-G-G-G, (7) 

z z z z z z 

with the first term the "direct" contribution and the sec- 
ond, the "exchange". We have left the z-dependence of 
the propagators explicit here, as this will be important 
for the calculation of NM effects. 



A. Self-consistent kernels 

The main idea behind this paper is to take the kernels 
obtained in LPT and replace in them the free system 
propagators with their effective counterparts that already 
include the effects of the reservoirs. The simplest SC 
kernel that one can come up with in this manner, which 
we denote S(")(z), is simply the SC equivalent of the 
sequential kernel of Eq. ©: 

S('^)(z) = G^G, (8) 

z 

in which the propagator "^" contains the effective system 
Liouvilhan W^^^^) = -Cg + T}-°''^{z) rather than just free 
Liouvillian £s ■ Expansion of this propagator shows that 
this self-energy reproduces the sequential term exactly, 
and provides an approximation to the direct cotunnel- 
ing term and all higher diagrams in which no contraction 
lines cross. This kernel is equivalent to a SC Born ap- 
proximation. 

The second SC kernel that we propose here reads 



I 1 I 1 1 I 

Y}^''\z) = G^G^rG^G^G^G, (9) 

z z z z 

where the propagator here contains Liouvillian 'W^^^ (z) = 
£s + S('')(z). This kernel explicitly includes not only a 
sequential-type diagram but also that of the fourth-order 
exchange term. Expansion shows that I](^'(z) contains 
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both exchange and direct cotunnehng contributions at 
fourth order — the exchange term exactly, and an ap- 
proximation to the direct term from the SCBorn part. 
At higher orders this kernel includes a far larger class 
of diagrams than does both with and without 

crossings. As we will see, this represents a significant ad- 
vance over the SCBorn kernel, S(")(z), not least because 
there is a large degree of cancellation between direct and 
exchange contributions. Lastly we also introduce the SC 
kernel 
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E('=)(z) = G-G + G-G^G^G 



-F 



-G — G ^ G ^ G, 



(10) 



which explicitly contains diagrams of both cotunneling 
types, as well as a mix of free and effective propagators. 
This kernel has the advantage that the sequential con- 
tribution is parameterised independently of higher order 
terms. Expansion shows that this kernel is exact up to 
fourth order. These three kernels all share the property 
that, upon expansion, they give a sum of diagrams with 
topologies the same as a subclass of diagrams of the exact 
kernel with no overcounting. 

Translated into an equation for the effective system 
Liouvillian, Eq. gives the following functional form 

W^'''>iz) ^Cs + f dt^iF('^)(>v(°)(0),wi), (11) 



where is some function (the form of which is ir- 

relevant for the current discussion, but can be found in 
Appendix [C|) , and where we have made explicit the inte- 
gration over bath frequency uJi. The corresponding func- 
tional form of Eq. ^ is the double integral 



with a similar equations for W^'^)(z). In all cases, we can 
make progress with these integral equations by introduc- 
ing the eigendecomposition of W on the righthandside 
(Appendix [B]). The integrals can then be evaluated ana- 
lytically as functions of the eigenvalues of W (Appendix 
[C|) and the resulting matrix equations can then be solved 
numerically by iteration. 



IV. COUNTING STATISTICS 

To calculate the current and its fluctuations we will 
employ the FCS formalism (26l . [271 as appropriate for 

master equation calculations ia[ifll2iii| 



The central 

object in this formalism is the counting-field-resolved Li- 
ouvillian W(x; z) = £s + 5^(xj z), in which processes that 
transfer n electrons to and from the leads are associated 
with counting factors e^*"-^° , where Xa is the counting 
field and a labels the lead in which counting takes place. 



Once in possession of the x-resolved Liouvillian, all zero- 
frequency cumulants of the current can, in principle, be 
calculated. Given the existence of this well-developed for- 
malism, the outstanding question for us to answer here is 
how to include the counting fields into our self-consistent 
kernels. 

Achieving this requires two modifications to definitions 
Eq. ©, Eq. ®, and Eq. Firstly, as in the LPT 

calculation of I, the bath contractions are modified to 
contain counting fields: 
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7fr(x) =7fre--«^^(^'-f^)>^-, (13) 



with Xcti J the counting field of lead ai, Sa = ±1, a factor 
given by the sign-convention for current flow in lead a; 
Pi and p2, Keldysh indices; and ~ ±1, an index in- 
dicating whether an electron enters or leaves the system 
at vertex 1 (see Appendix |^ and I for a full explana- 
tion of these symbols) . Secondly, the SC self-energies are 
not to be evaluated using propagator , but rather its 

z 

X-dependent counterpart " = " , which describes propa- 

x;z 

gation under the action of x-dependent kernel W(x; z). 
The equation for the x-dependent SC kernel of type (a) 
thus reads 



s^"^(x;^) 



G - Gg2(x), 



(14) 



where the propagator on the righthandside is evaluated 
with the full nonMarkovian x-dependent effective Liou- 
villian W^'^Hx; z)= Cs+ S('')(x; z) and where 



92 (x) 



(15) 



is the counting-field factor arising from the single bath 
contraction at sequential order. Similarly, for the other 
two schemes we have 



r 



-F 



SW(X;2) = G - Gg2(x) + G ^ G - G - G q^xix), 

x;z x;z x;z 



(16) 



and 



s(=)(x;z) = G-Gg2(x) + G-G - G - G 94d(x) 

z z x\z X\z 



-F 



(17) 



+G-G - G - Gg4x(x), 

z x;z x\z 

with counting-field factors 

g4D(x) = e'''°i«^^(fi-P^)'^°ie^''°2«2i(P2-P3)xo2^ (18) 

from the two bath contractions at fourth order. 

The justification for including the x-dependent self en- 
ergies S](x;2) on the righthandside propagators, instead 
of e.g. E(x = 0;z), comes from expanding these forms 
and comparing with an exact x-dependent LPT expan- 
sion. Moreover, including x in this fashion ensures charge 
conservation as discussed in the section IVl 



4 



Current correlations 



and the required second-order derivatives give 



Solving the self-consistent equations Eq. ([T^ , Eq. (IT51) , 
or Eq. pT)) for the full x-dependent E(z; x) as a function 
of z and x is most likely too demanding to be achieved in 
practice. Nevertheless, these equations may be used to 
generate the various quantities required to calculate cur- 
rent correlation functions. Here, we calculate the zero- 
frequency current correlation functions using the non- 
Markovian "current block" formalism of Refs. [lol[20l . [3l| . 
which does not require the full S(z; x), but rather its var- 
ious derivatives evaluated at z = x = 0. 

We work in a representation in which the elements of 
the system density matrix are arranged into a vector, 
such that Liouvillian >V(z) can be written as a matrix. 
The stationary state of the system, which we write as 
^stat _ 1^^^^ ^ggg Appendix iBjl . can be obtained as the 
right null-vector of the zero-frequency limit of W(z = 
0): W(O)lV'o)) = 0, and this we assume to be unique. 
Similarly, the left null- vector defined via = ((i/'o|W(0) = 
is the system "trace" vector [2§|. We further define the 
stationary state "expectation value" ((•)) = • iV'o)), 
and the projectors V — |7/'o))((^o| and Q = 1 —V. We 
then define 

J(X, e) = W(x, z = - ie) - W(x = 0, z = 0), (19) 

with the derivatives J' = c^x^lxe^O' ^ ^ '^eJ'l^^g^o' 
and analogously for higher-orders. We also require the 
zero-frequency pseudo-inverse 



n = hm St 



1 



£^0 ze -h >V(0, -ze) 



(20) 



The stationary average current can then simply be writ- 
ten as 

{I) = m), (21) 

with the zero-frequency shotnoise, S*, given by 

i^s = ij")) - 2ij'nj')) - 2(7) Uj')) - ij'nj)) 



(22) 

The final term in Eq. is the nonMarkovian correc- 
tion. Expressions for higher current cumulants can be 
found in Refs. [3l|, and obtained recursively but we 
focus here on the average current and shotnoise. 

Equations for the blocks required in the above expres- 
sions can be obtained by differentiating the definitions of 
the kernels Eq. ([13]), Eq. ([H]) or Eq. (H?]), and setting 
z, X — ^ 0. For the purpose of exposition, we consider just 
the (a)-type kernel here; the equations for the (6)-kernel 
are given in appendix |DJ The first derivatives of Eq. (|14p 
yield the equations 



J' 
j 



G 



G 



(il + j) 



—1 

G, 



(23) 
(24) 



iiJ")) = {{ G ^ G)) q'^ + 2{{G ^J'^ G)) q', (25) 

{{j')) = {{G ^ (zl + j) ^ G)) q', (26) 

In the latter case, we give expressions for the stationary 
expectation values of the blocks, rather than the blocks 
themselves, since these are simpler with many terms giv- 
ing zero in the expectation value thanks to the "leftmost- 
G rule" [131 ■ The evaluation of diagrams of the form 

G ^ M =i: G is discussed in Appendix [Cj Equations 
p3| and ((24)) each represent a set of linear equations for 
the matrix elements of J'' and respectively, which can 
easily be solved and the results substituted into Eqs.([25|) 
and (PS|) . 

Given the above formal developments, it is perhaps 
useful to give a brief overview of its practical applica- 
tion. Let us refer to the procedure for the (a)-type ker- 
nel. Firstly, all the basic elements of LPT are constructed 
as matrices, e.g. the free Liouvillian, the various tunnel 
operators, etc. Next the stationary kernel E(°)(0;0) is 
determined using Eq. ([5]). This equation is solved iter- 
atively with, for example, the standard LPT kernel as 
starting point. Once I]'°-'(0;0) is known, the station- 
ary state and pseudo-inverse TZ can easily be found. 
Equations ([23]) and ([24]) are then constructed, solved for 
J'' and j and the results being fed into Eqs. (P5|) and 
((26)) to obtain {{J")) and {{j')). These various elements 
are then combined in accordance with the expressions 
Eq. ((^T)) and Eq. ((^^ to give the average current and 
shotnoise. 



V. CHARGE CONSERVATION 

Before applying the method, we first demonstrate that 
it respects charge conservation. With charge conserved, a 
difference in instantaneous currents gives rise to a change 
in the accumulated system charge Q = Na\a){a\ with 
Na the number of electrons in state a (e = 1). In oper- 
ator terms we have Sq/q = —dQ/dt, where currents 
Sala are defined positive for electron fiow into reservoir 
a. Taking the stationary expectation value gives immedi- 
ately the relation Sq, (/„) = for the average currents. 
In terms of the jump super-operators of section dV Al this 
relation reads 



(27) 



with block J'^ is defined as in Eq. ((23)) , but with the lead 
index of the x-derivative made explicit. Similarly, for the 
shotnoises, charge conservation implies [11] 



(28) 



al3 



with the zero-frequency shotnoise correlator S'a^(O) be- 
tween leads a and (3. Using the lead-specific analogue of 
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Eq. (|22p . in terms of jump operators this relation reads 



(29) 



al3 



where the appHcation of Eq. (I?fl) makes the nonMarko- 
vian contributions cancel under the sum. 

We show that our approach obeys these relations by 
extending the rate-equation work of Ref. @ to the SC 
kernels studied here. Considering the E^") SC-sequential 
kernel as an example, explicit multiplication alows us to 
calculate various commutators involving system charge 
operator Q. Firstly, the commutator of Q with the full 
X-dependent kernel of Eq. reads 



>V(X),Q = Y.s^G^Gq'Jx) 



(30) 



(for brevity, we write here W(x) = W*°'(x; z = 0)). At 
X = 0, this relation yields 



(31) 



Differentiating with respect to ixp , setting X/3 ~^ and 



summing over /3 with lead- factors S/j we also obtain 
W(0),( 



S/sG ^ Jp ^ G,Q 



- IgTg q^piO) + G^J^G 

An identical calculation can be performed without the 
X-dependence of the propagator in W(x) which yields 



G - G g(x), 



= E*"^ - Gg:,(x). (33) 



This recovers Eq. pip and gives additionally 

r 



>V(0),gJ ,gj =^s„,s^G^Gg::^(0). (34) 
Comparison of Eq. ([32|) and Eq. ([34|) then yields 



Y,spG^J!,^G,Q 



Equation [21] allows us to write 



a a 

= 0, (35) 

where the first term vanishes due to the action of W(0) 
in the stationary expectation value, and the second due 
to the " leftmost- G rule" (H. Similarly, using Eq. ^ 
and Eq. ([35]) . the sum of shotnoises evaluates as 



X — ^ 

2_^Sj3G ^ Jp^G,i 



-2((|^[w,q] +^s„G- J-^-G^ QW-iQ |^[W,Q 



SaG ^ J'^ G 



(36) 



once again courtesy of the properties of W(0) and the 
leftmost- G rule. 



The mean currents and shotnoises therefore obey the 
properties demanded of them by charge conservation. 
An identical argument can be made for the other self- 
consistent kernels considered here. It is also relatively 
straightforward to extend this argument to higher cumu- 
lants and to other kernels. 



VI. RESULTS 



As an example application of this formalism, we con- 
sider a quantum dot with a single spinless level. The 
Hamiltonian reads H = ed^ d+H^es+V , with e the energy 
of the level when occupied and with lead and coupling 
Hamiltonians i?res and V as in Sec. [ill This model acts 
as an important benchmark for transport calculations, 
e.g. [3^ . as interactions are absent and exact solutions 
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eV/keT 



FIG. 1. Stationary current (/) through the single resonant 
level as a function of applied bias eV with system-lead cou- 
pling r_L = r_B — 5A;_bT. The main panel shows the current 
itself with results plotted for the three SC kernels discussed 
here, together with the exact solution and the fourth-order 
LPT results. For this coupling, the LPT solution breaks down 
when dot level is located around resonance with the chemical 
potentials. In contrast, all three SC results follow the exact 
solution closely. The inset shows the absolute value of the 
difference 51 between exact and each approximate solution. 
Further parameters were dot level position e = 20kBT, chem- 
ical potentials fiL = — ~ eV/2, and cut-off Xc = SOOfcsT. 



eV/ksT 




eV/ksT 

FIG. 2. Zero-frequency shotnoise S for the single resonant 
level as a function of bias. The main panel shows a compar- 
ison of the results from the self-consistent E'-''-' scheme with 
the exact results for three different couplings: Vr = Tl = 
1, 3, 5 ksT. The inset shows results for fourth-order LPT and 
the other two SC kernels for a coupling of Vh = Vl = bksT. 
Other parameters as in Fig. [T] 



are available from, e.g., scattering theory [13, H^- The 
exact results were summarised in I. 

In Fig. [T] we show the current calculated with our 
three SC kernels alongside the exact result and that from 
straightforward fourth-order LPT. Results are shown for 
a system-lead coupling of Tn = Tl = hksT. As is clear, 
at this coupling the fourth-order perturbative expansion 
LPT solution breaks down badly when the dot level is 
located around resonance with the chemical potentials. 
This is in line with the finding of Rcfs.[l3, [lO], which 
found the fourth-order ME calculation to be reliable 
across the full bias range for a coupling Tn/ksT < 1/2. 
In contrast, all three SC kernels give a very good account 
of the current. Scrutiny of the inset, which shows the dif- 
ference between exact an approximate solutions, shows 
that E(°) is the most accurate across the bias range, and 
that S*^"^ gives comparable results to the other kernels 
despite being much easier to calculate. 

Figure [2] shows results for the shotnoise of this model. 
Starting with the inset we see that that the LPT4 results 
once again break down around resonance. Moreover, we 
see that neither S'"' nor S^^^ provides a good description 
across the whole bias range. The results for E^'^) break 
down around resoncc in a fashion similar to the LPT 
solution and this can be attributed to the inclusion of 
unbroadened propagators in the kernel. Results at large 
and small bias for this kernel are nevertheless very good. 
Results for kernel E'"-' break down in a different man- 
ner; there is no abrupt behaviour around resonance, but 



the predicted value of the shotnoise is everywhere much 
higher than the exact solution. This difference is sig- 
nificant and remains even when the coupling is reduced. 
This error can be traced to the failure to includes di- 
agrams with crossings. In particular, at fourth order, 
direct and exchange contributions are of similar magni- 
tude, but with opposite sign. Failure to include diagrams 
with crossings therefore leads a large portion of the di- 
rect contribution remaining uncancelled, and this leads 
to the overestimated shotnoise. This problem is not man- 
ifest in the current because, for this model at least, the 
fourth-order exchange term docs not contribute anyway. 
Wc must therefore conclude that this SCBorn-typc kernel 
is unreliable for calculating shotnoise, and by extension 
higher counting statistics in this intermediate coupling 
regime. 

The main panel shows results for kernel E'^''-' and we 
see that this kernel does give a good account of the shot- 
noise across the bias range. For higher couplings, the SC 
results remain qualitatively correct except at low bias, 
when the system is in the cotunneling regime. Here the 
SC solution shows a spurious increase of the the shot- 
noise with decreasing bias. This arises from the conflict 
of requiring both that the sequential contribution vanish 
in this cotunneling limit, but also be finite to generate 
the correct direct cotunneling diagram under iteration. 
Clearly these two requirements can not be met at the 
same time, and this leads to a loss of accuracy in this 
regime. 
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VII. CONCLUSIONS 

In this paper we have described a method for cal- 
culating counting statistics using self-consistent kernels. 
Three different SC kernels were introduced and we have 
considered current and shotnoise for the single resonant 
level as an example application. For this model, all three 
kernels provide good results for the current, significantly 
better than could be achieved with straightforward per- 
turbative techniques. On the other hand, only S^''^ is 
able to provide good results for the shotnoise, in par- 
ticular when the dot level is around resonance with the 
chemical potentials. This is the most interesting regime 
since simple LPT performs well at both high bias and in 
the cotunneling regime. 

Although the above comparison is made for a non- 
interacting model, we expect the conclusions to hold 
broadly true for interacting models in the intermediate 
coupling regime considered here. It remains as future 
work to see how these SC techniques perform in strongly- 
coupled interacting systems in which genuine nonp ertur- 
bativc effects, such as the Kondo resonance [H,|33, come 
into play. 
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Appendix A: Diagrams 

Although we will leave much of the detail of the LPT 
expansion to I, it behoves us here to discuss those ele- 
ments relevant to the translation of our diagrams into the 
corresponding analytic LPT expressions. The (nonSC-) 
LPT sequential diagram of Eq. ^ reads 



S(2)(z) ^G-G, 



(Al) 



and this translates as 

I](2)(z) = (-l)Grf]s(^i)Gr7fr- (A2) 

Similarly, the fourth-order exchange term has the dia- 
gram 



,AVS I 1 I I 

E(^^'(z) = G-G-G-G, 



z z z 



(A3) 



which translates as 

e(4^)(z) = (-l)Gf 17s(z3)Gri7s(z2) 

xGrf^s(;.i)Gf7fr7fr- (A4) 

In these expressions G^ are Liouvillc-space tunnel ver- 
tices which add or remove electrons from the system, with 



multi-index subscript '1'= (^i,fci,ai) denoting direction 
^ = ±1, bath momentum ki and lead {a = L, R) indices, 
and with Keldysh index superscript {pi = ±1). These di- 
agrams include the free system propagator, which trans- 
lates as 

- ^ nsiz,n) = (A5) 

z Zm — L-S 

with free system Liouvillian Cs and with the subscript 
TO determined by the position of the propagator in 
the diagram. The frequencies Zm are given by z„j = 
z + 'Y^^=m+i ■'^i with Laplace variable z, n the order 
of the diagram (e.g. n = 2 for sequential order) and 
xi = —i£^i{uji -\- fj-ai), with reservoir frequency w/ and 
chemical potential /i^j . The reservoir contraction reads 

-/2r = %pi/(-eiPi^i) (A6) 

with Fermi function /, and index notation '1'= 
(— 'Ci, fci, ai)- Finally, the sign forefactors in Eq. ()A2p 
and Eq. (|A4[) come from a product of (— «)" for the nth 
order diagram and Wick sign (— l)^p with Np the per- 
mutation number of the digram (here Np = for the 
sequential term and Np ~ 1 for the exchange). In all 
diagrams, summations over all ^, a, and p, as well as 
integrations over all bath frequencies, are implied. As in 
the diagrams of [l^, one could label each tunnel vertex 
with appropriates indices, and each propagator with ap- 
propriate frequency index. However, since these are given 
uniquely by position in diagram, there is little advantage 
in doing so. We will, however, label the propagators with 
Laplace argument z, and/or x when appropriate, as our 
formal manipulations require explicit functions of these 
variables. 

The self-consistent equations of Eqs. ([8]), (|9]) and (fTO|) 
include diagrams as above but with propagator 

=^ V 

z Z,n - jOs - 2^{Z) 

where the self-energy ^(z) is the appropriate one for the 
calculation in hand, e.g. etc.. Moreover, for cal- 

culating the counting statistics, the relevant propagators 
are the x-dependent ones 

r^-^( (^^) 

2:„, - £s - S(z;x) 

In discussing x-dependent quantities, multiplying a dia- 
gram by counting- field factor q{x), corresponds to mul- 
tiplication of the corresponding super-operator by the 
same quantity (under the summation ofcoursc). For ex- 
ample, the full z and x dependent sequential self-energy 
diagram of Eq. ([Tl)) reads 



S('^)(z;x) = G - G<Z2(x), 



(A9) 



and this translates as 

(x; z) = Gf P'^^ ^'P'^';' Gf , 

where we have used the explicit form of 92 (x) from 
Eq. (USD. 
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1. Kernel expansions 



The eigendecomposition of the stationary kernel reads 



Insight into the approximations provided by our self- 
consistent kernels can be obtained by expanding the ef- 
fective propagators 



In this way, we can expand the SC self-energy T.^'^^z) of 
Eq. dH]) as 



G + G 



G 



— I 

G 



G 



(AlO) 



where we have made explicit the x variables in the 
propagators and, in doing so, have taken into account 
the (5-functions in the bath contractions. The expan- 
sion of therefore generates a diagram with the 
same topology as the direct cotunneling term. With x- 
arguments, the exact direct diagram reads 



G - G - 

z — xi z — xi- 



G 



G. (All) 



Comparison of Eq. (|A10p and Eq. (|Alip shows that the 
fourth-order term in the SC expansion provides an ap- 
proximation which differs only the x argument of the 
middle propagator. 



Appendix B: Matrix representation 

In this Liouville-space approach, superoperators 
such as the Liouvillian can be represented as finite- 
dimensional matrices and we often make use here of their 
eigendecomposition. Of particular importance is the de- 
composition of the X = 0, z = effective Liouvillian 
W — W(0,0+), from which the stationary properties of 
the system may be calculated. Since W is nonHermitian, 
it has non-adjoint left, {{ipal, and right, \ipa)), eigenvectors 
defined by the equations 



W\2Pa))=Wa\%)) 



(Bl) 
1 with 



with eigenvalues Wa and index a — 0, 1, 
N the dimension of W. Taken together the left and right 
eigenvectors form a bi-orthonormal set: {{ipa\4'a')) = ^a,a' 
and we have the closure relation in Liouvillian space 



(B2) 



^ Wa|V'a))((V'a| 



(B3) 



Assuming that the system has a unique stationary state, 
one (and only one) of the eigenvalues will be zero, wq = 0, 
and the stationary state itself is given by the correspond- 
ing right eigenvector pf'^^ = \iPq)). The dual vector {{ipo\ 
is then the "trace vector" with elements unity at posi- 
tions mapping to populations and zero at those mapping 
to coherences. The remaining eigenvalues of W have non- 
positive real parts and, if complex, occur in complex con- 
jugate pairs. 



Appendix C: Integrals 

We now discuss the analytic evaluation of the integrals 
required here. For convenience we set ksT = 1. The 
z = integral equation for S(°^(z = 0) reads 



-Pifi-^iPi^i) 



o+ + 26(wi + Moi)-w('^Ho) 



Gf . 



Using the eigendecomposition for W^") (0), as in Eq. (jB3|) . 
we write this as 



e(")(o) 



p,Gf\^a))Ua\G^,' 

/(-^IPlwl) 



duji- 



= -2^PiGf |^,))((^,|Gf 4f{(l^io., 



iWa), 



(2) 

which defines the sequential integral, Ip (A). This 
integral can be evaluated as in I with the intro- 
duction of a Lorentzian cut-off function 'D(ll)) = 

Xl,/{{u;-u;or 
convenience, as 



} , with band-centre ujq chosen for 



duj- 



_.l;,pA) + |*(A)^ 



Here 0(a;) 
with g{uj) 



5(w) 



g{-uj) - g{uj + iXc) - g{-uj + iX^)) , 
I and \I' the digamma function. 



As discussed in I, evaluation of fourth-order terms requires the integral 

.f{Pl(^l)f{P20J2) 



■ {iO+ + uji+uj2~ A2)(iO+ + cji - A3) ' 



(CI) 



In straightforward Liouvillian Perturbation theory, the Xi are real; for self-consistent calculations, they can be complex 
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and this means that we need the result for this integral in a slightly more general form than given in I. We have 

P'PH^2,h)=PiP2G{X2,X3)+piG{X3)+piP2H{X2,X3)+ const (C2) 
with "const" an unimportant constant. The first two part functions are 

G(A2, A3) = 27r2(i + 6(A2)) {/i'i^(A3) - /i'i'(A32)} + .g(-A3)(/.(A32) 

-55(-A3 + iXcmX32 - iXc) - i<?(A3 + iXcmX32 + iXc), (C3) 

and 

G(A3) = -7r2/|2)(A3), (C4) 

with Ay = Xi — Xj and Bose-Einstein distribution function b{x) ~ (e^ — 1) ^. The third contribution is more 
complicated and reads 

°° r 1 1 2 1 

HiX2,X3)^ 9{-X2+ink)\- .y ^ . , + ■ ^ .y ^ . , - , ^ . , 

^^^g L A32 - iXc + iTTk A32 + iXc + tnk A32 + ink J 



,3 

— g{—^2 + iXc + ink) 
+ !]{^2 + ink) 



1 1 2 



A32 + ink A32 + 2iXc + iTrfc A32 + iXc + iTr/c 
1 1 2 



-A32 — iXc + ink — A32 + iXc + ink — A32 + ivrfc 

f 1 1 2 

fl(A2 + + iTrfc) <^ h L (C5) 

^ \ -A32 + i^A: -A32 + 2^Xc + z^A: -A32 + + ivrfc ^ ^ ^ 



I 1 

It remains to discuss the translation and integration of diagrams of the form G ^ M ^ G, such as occur in Eq. 
for example. The analytic expression of this diagram reads 

Introduction of the eigendecomposition for the stationary kernels and the employ of a partial fraction decomposition 
allows us to write this as 

-p,Gf\i^a))mM\MU,,\Gl^ ( du J^-^'P'^^'U \ ^ \ ^ 

J Wa - Wa' \0+ + iCl (l^l + fJ-ai ) - Wa 0+ + 1^1 (^1 + Pai ) - Wa' 

Identifying the sequential integrals, this becomes 

- 2np,G\^\i:,))i^a\MMi^,,\Gf 4f (/^^i + " 4f (Mo, + ^^/;aO _ (^7) 



In the case that Wa = Wa', the differential quotient is recognised and the integral part replaced with d/ dwalp^ {p.ai 
iwa)- The evaluation of similar cotunneling-type terms follows analogously. 

Appendix D: Jump-operator equations for the E'*' SC kernel 

The four equations corresponding to Eqs p3ll26[) for the current blocks of E'^^ are: 



I 1 I 1 I I I 1 I 1 I I 

J' = G^Gq'2 + G^G^G^Gq'ix + G^J'^G + G^J'^G=^G^G 



I 1 I 1 I 1 

G^G^J'^G^G + G^G^G^J'^G (Dl) 



J^G^[il + J]^G + G^[il + J]^G^G^G 



G^G^(il + j^^G^ + G^G^G^(il + j^^G (D2) 
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I — I I — — I 



I 1 1 ^ 1 I 1 I 1 1 ^ 1 

2{{G^J'^Gq'2 + G^J'^G^G^Gq'4x + G^G^J'^G^Gq'4x + G'^G^G^J'^G q'^x)) 

(D3) 



(( j')) = ((G - (^l + j) - G + G - (i.l + ^ G ^ G ^ G q'^^ 



1 I — — : — ^ — I 



+ ((G = G = (il + j) = G - G + G - G - G - (tl + j) - G g^^)). (D4) 
The corresponding blocks for E^"^) follow analogously. 
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